Rare Events, Adaptive Umbrella Sampling, and Metadynamics

This notebook tries to illustrate the concept of a rare event, the logic behind employing an umbrella potential, and Metadynamics (both the plain and the well-tempered version). Everything is implemented trying to keep things as simple as possible. The notebook relies on numpy, matplotlib, scipy, and tqdm.

Everything will be performed on a small toy system. The system we consider is 1D. We use the following functional form for the potential energy:

\begin{equation} U(x) = ax^4 - 6x^2 + x \end{equation}

with $a=0.5$, $b=-6$, and $c=1$. It follows that we have a very simple expression for the force:

\begin{equation} f(x) = -\frac{dU(x)}{dx} = -2x^3 + 12x - 1 \end{equation}

The potential energy has the following look:

Then we define a simple System class. It's nothing more than a container with the initial coordinates and an attribute traj that will store the trajectory

We can now initialize our system. We initialize it in the basin on the right (x0 = 2.4), with no initial velocity (p0 = 0.0), with mass equal to one (m = 1), and with a small $k_bT$ (kbT = 1e-2)

We can take a look at it

Here we define a simple class that represents a molecular dynamics engine. We simulate our system in the NVT ensemble, by means of the Langevin equations of motions (Langevin dynamics). We integrate the equations of motions using the BAOAB method.

\begin{align} p(t+\frac{1}{2}\delta t) &= p(t) + \frac{1}{2}\delta t f(f) &B\\ x(t + \frac{1}{2}\delta t) &= x(t) + \frac{1}{2}\delta t p(t + \frac{1}{2}\delta t)/m &A \\ p'(t+\frac{1}{2}\delta t) &= \exp(-\gamma \delta t)p(t + \frac{1}{2}\delta t) + \sqrt{1-\exp(-2\gamma\delta t)}\sqrt{mk_bT}G &O \\ x(t+\delta t) &= x(t + \frac{1}{2}\delta t) + \frac{1}{2}\delta t p'(t + \frac{1}{2}\delta t)/m &A \\ p(t + \delta t) &= p'(t+\frac{1}{2}\delta t) + \frac{1}{2}\delta t f(f + \delta t) &B \end{align}

Here $G$ is a Gaussian random variable with zero mean and unit variance, and $\gamma$ is the frictional coefficient.

The class is initialized with some MD settings like the gamma constant of the Langevin dynamics, the time step, the number of steps to do, and the function that evaluates the force. It has a method run that runs the simulation and stores the coordinates in the given System.

High barrier

Now we can initialize the system and the engine. We first consider the case of a high barrier. "High" here refers to the size of the barrier compared with $k_bT$, i.e., we set $k_bT$ to a small value.

If we plot the trajectory, we see that the system visits only the basin on the right

We can also visualize the trajectory directly. This may take some time, and more so if the trajectory gets longer.

We are clearly stuck inside one basin. The probability distribution of our collective variable (the only variable present in this system, to be fair) reflects precisely this behaviour

Of course the corresponding free energy surface cannot be the right one

Small barrier

Now let's consider the case of a small barrier, i.e., a barrier that is not much bigger than $k_bT$. We then set $k_bT$ to a bigger value

Now we can cross the barrier

Still, the probability distribution is wrong. This is related to the sampling time. With "infinite" sampling, the probability distribution in this case of a small barrier would be a good estimate of the real one. With limited sampling time, we end up with a bad probability distribution. The corresponding free energy surface will also be wrong, of course.

More sampling time (small barrier)

We still keep the barrier low, but we improve the sampling by a factor of 10. In this way, we should be able to observe a large number of transitions, and obtain a good estimate of the free energy. This situation reflects the ideal case of "infinite" sampling.

Here, with a lot of simulation time, we actually see a lot of transitions, and the system spends more time in the basin on the left, as it should do.

Umbrella Sampling

Let us consider again the case of a small barrier. We want to find a way to cross the barrier. The idea of Umbrella Sampling is to put an "umbrella" potential that helps the system to cross the barrier (note that this is not the standard way of doing Umbrella Sampling in real cases, where harmonic potentials are used to restrict the exploration to specific intervals along the CV). Let us try a guess for the umbrella potential. Given the shape of the potential energy, an umbrella with the shape of a bell should work. Not again that this kind of reasoning can be made only given that this system is very simple.

This seems enough to cross the barrier. Let's add the corresponding term to the force

We run the simulation again with this modified potential

In this way we can actually cross the barrier

Now that we have crossed one basin, we can do the trick again. Let us guess the next gaussian potential.

With this potential, we should achieve an almost free diffusion in the available space, meaning that we are sampling everything correctly.

Now we can diffuse (almost) freely between the two basins. They are both completely accessible.

Note that the probability distribution of $s$ is not uniform over the whole interval, as the bias potential we are applied does not perfectly makes the resulting potential energy "flat" in the region of interest. Therefore, we have peaks in the probability distribution corresponding to the small bowls in the potential energy. Note also that this is a probability density of a biased system.

As this is the probability distribution of a biased system, we can't recover the free energy of our true system by simply taking the (minus) logarithm of this quantity. Let's try to do that

Instead, we must use the main result of Umbrella Sampling

\begin{equation} F(s) = F_B(s) - B(s) + C_B \end{equation}

where we don't care about $C_B$. If we do that we get a good free energy. By sampling more and more our estimate of the free energy would become even better.

In this case we have guessed the umbrella potential for each of the two basins. An alternative is to use Adaptive Umbrella Sampling, where instead of guessing the potential one runs a small MD simulation, computes the $P(s)$ in that simulation, converts it to a $F(s)$, and uses $B(s)=-F(s)$ as a bias potential for another small simulation. Then repeats the process, iteratively building up $B(s)$. We do not implement such a procedure here. Instead, we skip directly to metadynamics.

Metadynamics

Let us now use metadynamics (not the well tempered version). In order to do so we must modify our MDEngine class and implement a simple version of metadynamics.

Now we look at the metad.history in order to track every added Gaussian, and we plot them over the potential energy. We see that the bias has covered all of our potential energy, so that we have free diffusion of the system, after an initial transient

Let's see the diffusion

Finally, we obtain the free energy surface directly from the bias

Now for the well-tempered version. We can activate well-tempered metadynamics by providing our Metadynamics class with the bias factor and the temperature of the system. The algorithm will add Gaussians, but the height of each Gaussian is rescaled according to the bias that has been already deposited. This should provide us with a smoother filling of the free energy. Furthermore, the Gaussians will eventually have "zero" height, so that, after a transient time, the system will not explore the free energy further. If the biasfactor is too low, then, there may be barriers that can not be overcome.

We can look again at the diffusion in CV space and the free energy